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From an extensive calculation of static properties of a trapped Fermi superfiuid at zero tempera- 
ture using a density-functional formulation, we demonstrate a universal behavior of its observables, 
such as energy, chemical potential, radius etc., over the crossover from the BCS limit to unitarity 
leading to scaling over many orders of magnitude in fermion number. This scaling allows to predict 
the static properties of the system, with a large number (~ 10 5 ) of fermions, over the crossover with 
an error of 1-2%, from the knowledge of those for a small number (~ 10) of fermions. 



PACS numbers: 03.75.Ss, 71.10.Ay, 67.85. Be, 05.30.Fk 



I. INTRODUCTION 

The Pauli principle leads to an effective repulsion in 
identical fermions which could dominate the physical 
interaction and control the dynamics. The dominance 
of Pauli repulsion is responsible for the stability of our 
hadronic universe. When this happens the system ex- 
hibits universal behavior practically independent of or 
weakly dependent on the physical interaction. One clas- 
sic example [l[ of this is found in the s-wave quartet 
nucleon-deuteron system with three spin-parallel nucle- 
ons controlled by Pauli repulsion. The calculated scat- 
tering length for different nucleon-nucleon interactions is 
essentially 6.3 fm, whereas the doublet s-wave scattering 
length (not controlled by Pauli repulsion) for the same 
interactions varies from — 2 fm to 3 fm pj, This uni- 
versality is prominent in the limit of zero Fermi-Fermi 
interaction in the Bardeen-Cooper-Schrieffer (BCS) the- 
ory of superfiuid fermions leading to universal properties 
of low-temperature superconductors [H , of cold neutron 
matter and neutron star 0, Q , and of a trapped Fermi 
superfiuid kj at zero temperature. This universality also 
manifests 0, d, S EH in a trapped Fermi superfiuid at 
unitarity as the Fermi-Fermi scattering length a goes to 
infinity (a — > — oo). In this limit, though the physical in- 
teraction is nonzero, the only interaction scale (scattering 
length) disappears and the system acquires universal be- 
havior. This interaction scale is also absent in the BCS 
limit (a = —0) with universal properties. Similar univer- 
sality is found in other fermionic systems with large \a\ 

The crossover from weak-coupling BCS limit to uni- 
tarity [131 has been a very active area of research 
@, 0, Ild.ll3j after the experimental realization [3, EH of 
this crossover in a trapped dilute Fermi superfiuid near 
a Feshbach resonance. Using a complete numerical sim- 
ulation of a density- functional (DF) formulation we 
show that the deviation from universality of a trapped 
dilute Fermi superfiuid over the crossover is orderly and 
the system continues to exhibit nearly universal behav- 
ior possessing useful scaling relations over many orders 
of magnitude involving energy, chemical potential, radius 
and number of Fermi atoms (N). In this crossover region 
the Pauli repulsion dominates over the physical interac- 



tion leading to the universal behavior. 

In Sec. II we present the density functional formula- 
tion that we use in this study. In Sec. Ill we present 
the numerical results and establish the universal behav- 
ior of a trapped Fermi superfiuid in the BCS-unitarity 
crossover. Finally, in Sec. IV we give some concluding 
remarks. 



II. DENSITY FUNCTIONAL FORMULATION 

To study the universality, we use a Galilei-invariant 
DF formulation for the crossover of a trapped two- 
component Fermi superfiuid [17] , equivalent to a hydro- 
dynamical model with the correct phase- velocity relation 
p v = KS70/{2m), where v is the superfiuid velocity, m 
the Fermi mass, and the phase of the order parameter 
^(r) at position r, which satisfies (Eq. (35) of [l7|, but 
with a distinct g(x) consistent with the known small- a: 
behavior of energy of a uniform Fermi gas) 
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V 2 + U + fi(n, a) 



^{n,a) = — (37r 2 n) 2 / 3 5 (n 1/3 a), 
2m 



*(r)= W) *(r), (1) 
(2) 



fJ-o 



dr 
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jV<F| 2 -I- U^ 2 + ^(n,a)* 2 



where J $> 2 (r)dv = N, and 



g( x ) 



(l-Pix + /3 2 x 2 )' 



(3) 



(4) 



with xi = M(3^ 2 ) 2 / 3 , X2 = 300, fa = 40, 2 = X2 /(l - 
C), fi(n,a) the bulk chemical potential, /zrj the chemical 
potential for the trapped system, n — "F 2 the density 
of atoms, and U — muj 2 r 2 /2 the harmonic trap of fre- 
quency lo. (Different parametrizations of g(x) were used 
in Ref. |l8|.) Here we take C = 0.44 consistent with 
Monte Carlo calculations 0, [l3[ and experiments of 
a uniform Fermi superfiuid at unitarity. The parame- 
ters Xi,X2,/3i and 02 are chosen so that the model (i) 
agrees with the fixed-node Monte Carlo (FNMC) 0] and 
Green-function Monte Carlo (GFMC) 0] results for the 
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energy of a trapped superfluid at unitarity (for N < 30) 
and over the crossover \A [l7| . (ii) provides a smooth inter- 
polation between the energies of a superfluid at the BCS 
and unitarity limits (TtI |21| , and (iii) satisfies the known 
BCS limit [22j of the bulk chemicaj. potential fj,(n, a) (two 
lowest-order terms of Eq. (1) of [23j). 

In Ref. [l7| we used the following simple expresion for 
the function g(x) in place of that given by Eq. (|4"|): 



1 



I -fix' 



(5) 



with x = 207r/(37r 2 ) 2 / 3 , and (3 = x /(l - 0- Tnis choice 
does not satisfy the known weak interaction BCS limit 
\im x ^ g(x) -> 1 + 47ra/(37r 2 ) 2 / 3 [13, El, HI- Equation 
([5]) was used in Ref. QjJ as a simple model to satisfy the 
BCS-unitarity crossover Monte Carlo results [2l[ for en- 
ergies of trapped fermions. The bulk energy (and chem- 
ical potential) of a Fermi gas in the weak-couping BCS 
limit can be written as an expansion in the parameter 
x = n 1 / 3 a as discussed in Ref. [23|. Lee and Yang 
|22| calculated the coefficients of the series in a specific 
model of interacting fermions as quoted in Eq. (1) of 
(23)]. These coefficients can be related to the coefficients 
in an expansion of g(x) in the small x limit. The x- 
independent constant in g{x) and the coefficient of the x 
term \ = 47r/(37r 2 ) 2 / 3 in this series are reasonably model 
independent. However, coefficients of x 2 will be depen- 
dent on the interaction model. In Ref. [TtJ we changed 
the coefficient of the x-term from that proposed in [23| , to 
simulate the effect of the (unknown) higher order terms 
in the x-expansion of g(x). The coefficient of the x-term 
can be kept at the Lee- Yang value 47r/(37r 2 ) 2 / 3 provided 
we include a x 2 term in the expansion of g(x) and we 
have done this in Eq. However, the coefficient of 

the present x 2 term does not and should not agree with 
that of the known expansion quoted in Eq. (1) of [23l |. 
This is because this coefficient now simulates the con- 
tribution of higher order terms. Although the analytic 
expression g(x) is modified here from that given by Eq. 
© to (fj| , from a numerical (calculational) point of view 
the change is negligible as we can see from Fig. [TJ where 
we plot the two parametrizations of the function g{x) 
given by Eqs. ((4]) and ([5]). In Ref. [l7[, we used the 
function g{x) of Eq. ([5]) to reproduce some Monte-Carlo 
results [2l| for energy of a trapped-fermion system over 
the BCS to unitarity crossover. We have now performed 
the same calculations with the function g(x) of Eq. ((4]) 
— illustrated below in Fig. Q] (b) — and we verified that 
the results for energy remain practically unchanged. 

The gradient term in Eq. ([T]) provides a correction to 
the local density approximation (LDA) [24[ obtained by 
setting the gradient term to zero. LDA is a good approx- 
imation for a large N, when the bulk chemical potential 
jj,(n, a), a positive term responsible for Pauli repulsion in 
the system even for attractive (negative) a, is very large. 
The gradient term is consistent with the hydrodynamic 
flow of paired fermions of mass 2m @, [TtJ ■ To study the 
scaling of the solution of Eq. (fl]) we note that in the BCS 
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FIG. 1: (Color online) The function g(x) vs. x as defined by 
Eqs. © and ©. 



and unitarity limits the bulk chemical potential Lt(n, a) 
has, respectively, the following simple forms [j| H, 
h 2 {3n 2 n) 2 / 3 /{2m) and C?i 2 (37r 2 n) 2/3 /(2r7i), which shows 
the scaling n(n, a) ~ n 2 / 3 in both limits. In these limits 
the energy functional is given by Q 



E 



dr 



8m 1 



U^ 2 



3ft 2 £ 
20m 



(3^2)2/3^10/3 



(Ey) + (Spot) + (^Fermi), 



(,6) 
(7) 



where £ = 1 in the BCS limit and £ = ( = 0.44 at 
unitarity [1, HI], and (E v ), (E pot ), and (E Fcrm {) are, re- 
spectively, the expectation values of the three terms in 
Eq. ^ . This "analytic" dependence of E on ^ leads to a 
simple virial theorem, postulated and studied experimen- 
tally in connecting (Ey), (Spot), and (SWim) at the 
BCS and unitarity limits. For the exact ^S>, energy E is 
a minimum. To derive the virial theorem, we take the 
Cartesian system r = (x, y, z), and the norm-preserving 
scaling transformation ^(x,y, z) — > y\ij){Xx,y,z). The 
condition of minimum energy is \dE / d\\\ = \ = 0. The A- 
dependent part of energy now becomes E\ — \ 2 {E^) + 
A~ 2 (Sp Qt ) + A 2 / 3 (Spcrmi) , where the suffix x denotes x 
component. The minimization condition yields (£p 0t ) = 
(Sy) + (SFermi)/3- By summing over three components, 
we have the virial theorem (E pot ) = (Sy) + (Spermi), 
or equivalently, E — 2(E pot ) at the BCS and unitarity 
limits. The deviation in percentage from the virial theo- 
rem in the crossover can be estimated by the percentage 
defect function t] = 100 (E - 2(E pot ))/E. 

In the LDA one has the following analytic solu- 
tions for the energy, chemical potential, and mean 
square radius of the trapped system [24j: E/(Kv) = 
(37V) 4 / 3 V?/4,/io/(M - (3A0 4 / 3 V?A (r 2 )/(h/mu;) = 
(3) 4 / 3 A rl / 3 v / C/4. The absence of the parameter a in 
n(n, a) of Eq. |(5J) in the BCS and unitarity lim- 
its leads to the following properties: (i) the scaling 
fi(n, a) ~ n 2 / 3 for the uniform Fermi superfluid, (ii) scal- 
ing E / (hu)) , ijLq / (tuv) ~ N 4 / 3 and (iii) the virial theorem, 
E = 2(E pot ) for the trapped Fermi superfluid. These 



3 



universal properties lead to predictability of the Fermi 
superfluid in the BCS and unitarity limits, e.g., predict- 
ing the energy for a large number of fermions from the 
knowledge of that for a small number of fermions ob- 
tained by accurate "exact" calculation. This is important 
as exact calculations for large systems are very difficult, 
if not impossible. 

To extend the above predictability over the crossover, 
where analytical (LDA) results are not available, we 
study the deviation from the universal properties of a 
Fermi superfluid in the crossover region. We find that 
the deviation is orderly, which allows us to restore a uni- 
versal and predictable behavior of the trapped Fermi su- 
perfluid, so that one can predict the properties of a large 
system from a knowledge of those for a small system. 



III. NUMERICAL RESULTS 

We solve Eq. (fTJ) by transforming it to time-dependent 
form by replacing /i by a time derivative. We express 
energy variables in units of hu>, length in \fhj (muj) and 
time in w . The resultant equation is then discretized 
by the semi-implicit Crank-Nicholson algorithm using a 
typical space step 0.04 and time step 0.001 and then 
solved by imaginary time propagation. The chemical po- 
tential fio is then calculated via Eq. ([3]), energy in the 
BCS and unitarity limits via Eq. ([6|), and the energy in 
the crossover region via a numerically constructed energy 
functional. In Table [J we display energy per particle in 
units of Ep = (SNy^hu; (the Fermi energy of an ideal 
Fermi gas at the trap center) or E/[huj(3N) 4 / 3 ] along the 
crossover for different N and a. In Table HT1 we report the 
respective chemical potentials fj,o/[hu;(3N) 4 / 3 ]. 

To understand the universal nature of E and /io of 
Tables U and [TT] for different N, we study the nonlin- 
ear input /z(n,a) to Eq. (p}. In Fig. [2] (a) we plot 
4ir 2 nn(n 7 a)/(h 2 n 2 / 3 ) vs. lain 1 / 3 for a < 0. In the BCS 
(jajn 1 / 3 — * 0) and unitarity (lain 1 / 3 — * oo) limits one has 
the perfect scaling, /i(n,a) ~ n 2 / 3 , with deviation from 
this behavior in the crossover region. Next, to study 
the behavior of E and /io of the trapped Fermi super- 
fluid, we plot in Fig. [5] (b) the numerically calculated 
E/N 4 / 3 and fJ, /N 4 / 3 vs. N for a = -0, -0.5 and -oo 
and the FNMC @, S and GFMC 0] results. Perfect 
scalings, /j.q,E ~ N 4 / 3 , are observed for large N in the 
BCS and unitarity limits, with deviation in the crossover 
region and for small N. A careful analysis of the fiQ,E 
data of Fig. [2] (b) for small N reveals in the crossover 
region the average scalings {E/N 2 / 3 — 0.37) ~ N 2 / 3 
and (no/N 2 ? 3 - 0.27) - TV 2 / 3 . If these scalings were 
really perfect then plots of (E/N 2 / 3 - 0.37)/iV 2 / 3 or 
(Vo/N 2 / 3 - 0.27)/7V 2 / 3 vs. |a| for a fixed N would lead 
to universal curves independent of N. 

To quantify the deviation of the Fermi superfluid in the 
crossover region from universality, we study the deviation 
of our result for the trapped Fermi superfluid from the 
virial theorem, E = 2(E pot ). For this purpose we plot 
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FIG. 2: (Color online) (a)Dimensionless bulk chemical poten- 
tiaf 47r 2 mu(n, a)/(ft. 2 n 2 / 3 ) of a uniform Fermi gas vs. lain 1 / 3 
for a < 0. (b) Chemical potentiai and energy iio(N,a)/N 4 / 3 
and E(N,a)/N 4 / 3 of a trapped Fermi gas vs. N for a = 
-0, -0.5, -oo. The energies of FNMC 0,1 and GFMC [Ij 
calculations at unitarity are also shown. (Energies and scat- 
tering lengths are expressed in oscillator units.) 
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FIG. 3: (Color online) (a) The plot of 77 = W0[E(N,a) - 
2{E pot )]/E(N,a) for different N vs. |a|; and (b) n = 
W0\E(N,a) - 2{E pot )]/E{N,a) vs. \a\, a = N 1/6 a. Scat- 
tering lengths are expressed in oscillator units. 
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TABLE I: Dimensionless energies E/[hu(3N) 4/3 ] of a trapped 
Fermi superfluid along the crossover. The last two columns, 
{4} and {10}, give the predicted energies for N = 10 em- 
ploying scaling (|8j) using the energies for N' = 4 and 10, 
respectively. 



a\N = 


10 


10 2 


10 4 


10 5 


10 5 {4} 


10 5 {10} 


-0.001 


0.2656 


0.2539 


0.2501 


0.2499 


0.2428 


0.2470 


-0.01 


0.2652 


0.2531 


0.2478 


0.2459 


0.2393 


0.2434 


-0.1 


0.2550 


0.2364 


0.2146 


0.2046 


0.2027 


0.2049 


-1 


0.2079 


0.1869 


0.1740 


0.1714 


0.1731 


0.1740 


-10 


0.1897 


0.1731 


0.1670 


0.1665 


0.1686 


0.1693 


-100 


0.1875 


0.1714 


0.1662 


0.1659 


0.1682 


0.1688 



TABLE II: Same as Table U for the chemical potential /io. 



a\N = 


10 


10 2 


10 4 


10 5 


10 5 {4} 


10 5 {10} 


-0.001 


0.3449 


0.3361 


0.3333 


0.3331 


0.3284 


0.3312 


-0.01 


0.3443 


0.3351 


0.3304 


0.3281 


0.3237 


0.3263 


-0.1 


0.3310 


0.3139 


0.2880 


0.2749 


0.2743 


0.2753 


-1 


0.2670 


0.2474 


0.2326 


0.2290 


0.2315 


0.2314 


-10 


0.2409 


0.2277 


0.2225 


0.2219 


0.2249 


0.2246 


-100 


0.2376 


0.2254 


0.2215 


0.2212 


0.2242 


0.2241 



the function 77 = 100(£ - 2(E pot ))/E vs. \a\ for different 
N in Fig. [3] (a). The shift of the curves with N clearly 
shows that the results are dependent on a and N. We 
find that the maximum deviation from the virial theo- 
rem in the crossover region is quite small and is about 
4%. A careful examination of Fig. [3] (a) shows that this 
dependence on \a\ is linear in log scale and can be in- 
cluded by plotting 77 vs. \a\,a = N 1 ' 6 ^ as can be seen 
from Fig. [3] (b). In Fig. [3] (b) we find that all curves for 
77 have collapsed essentially on a single curve. However, 
the curve for TV = 10 is a bit different from others, which 
confirms that there is a deviation from universality for 
smaller N. This \a\ — |A rl / 6 a| dependence of properties 
of a trapped Fermi superfluid is quite universal as we see 
in the following. 

Now we study the universality in E and /io of the 
trapped Fermi superfluid obtained from a solution of Eq. 
((I]). We plot the function 5 = [E(N, a)/N 2 / 3 -0.37]/N 2 / 3 
vs. \a\ for different N in Fig. [4] (a) [We recall that the 
universal nature of the function 5 was obtained from an 
analysis of results in Fig. [2] (b)]. We find that as in Fig. 
[3] (a) the results for different N are distinct. If the sys- 
tem were really dominated by universality, a plot of the 
scaled quantities k = [E(N, a)/N 2 / 3 - 0.37}/N 2 / 3 vs. |a| 
would lead to universal curves. This is indeed found in 
Fig. g] (b), where we also included the FNMC results [7] 
for small N = 4, 8. For small N there is some deviation 
from universality which disappears for N > 10. We note 
that in Fig. H(b) the N = 30 and in Fig. H](b) the large- 
N FNMC and GFMC data lie on the universal curves. 
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FIG. 4: (Color online) (a) S = [E(N, a)/N 2/3 - 0.37]/A 2/3 vs 
|o| for different N for a trapped Fermi gas from Eq. ([lj. (b) 
K = [E(N,a)/N 2/3 - 0.37]/iV 2/3 vs |a| for different JV from 
Eq. CEJ and FNMC method HH- (c) v = [fi (N,a)/N 2/3 - 
0.27]/iV 2//3 vs \a\ for different N. (Energies and scattering 
lengths are expressed in oscillator units.) 



These universal curves yield the simple formula for the 
energy of a large system with atoms in terms of that 
of a small system with N' atoms (N >> N') 



E(N, a) = 



4/3 E(N',a)/N' 2 / 3 -^ 

AT'2/3 



(8) 



with 7 = 0.37, and where we have neglected the 
small constant 0.37 compared to the large quantity 
E(N,a)/N 2 / 3 . When N and N' are both large, formula 
© becomes E(N, a) = E(N',a)(N/N') A / 3 . We note 
that, in Eq. ([5]), energies are to be considered for the 
same a and not a. We indeed calculated the energies 
for A^ = 10 5 atoms using the data for N' — 4 and 10. 
The predictions so obtained for N' — 4, 10, listed Table 
HI compare well with the calculated results within about 
3% and 2% errors, respectively. In Fig. [4] (c) we plotted 
scaled chemical potentials vs. \a\ for N = 10 to 10 6 and 
find that they all lie on the same universal curve. In this 
case also a prediction of chemical potential for large N 
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FIG. 5: (Color online) (r) /N 1 ^ 6 of a trapped Fermi superfluid 
vs. TV for different a — — 0, —0.5, — oo. (Lengths are expressed 
in oscillator units.) 



using the same for a small AT' can be made through the 
scaling formula {8]) but now with 7 = 0.27 as obtained 
from an analysis of Fig. [2] (b). The predicted jUo for 
AT = 10 5 using the data for AT' = 4 and 10, listed in 
Table [Til compare well with the calculated results within 
an error of less than 1.5%. 

Because of the universal behavior due to the domi- 
nance of Pauli repulsion, other observables of a Fermi 
superfluid are also correlated with energy and atom num- 
ber. For example, the root mean square radius (r) has 
the scaling (r) ~ N 1 / 6 in BCS and unitarity limits. In 
Fig. \E\ we plot {r)/N x / 6 vs. A" for different a. Again 
there is a slight violation of this scaling for small AT and 
in the crossover region, which could be remedied by a fine 
tuning of the type Nevertheless, this clearly shows 
that the radius of a large trapped Fermi superfluid is pre- 
dictable from the knowledge of radius of a small system 
not only in the BCS and unitarity limits but also along 



the crossover. In the case of nuclei, where Pauli repulsion 
plays an important role, correlations exist among binding 
energy, radius and the number of nucleons [26| . 



IV. CONCLUSION 

In conclusion, from an extensive numerical study of 
the static properties of a trapped two-component Fermi 
superfluid using a Galilei- invariant DF formulation [l7T |. 
equivalent to a generalized hydrodynamic formulation 
with the correct phase- velocity relation [f|, we establish 
that, because of the dominance of the Pauli repulsion, 
the trapped Fermi superfluid has a universal behavior 
not only in the BCS and unitarity limits but also in the 
crossover region. This allows a prediction of the static 
properties (energy, chemical potential, rms radius etc.) 
of a large Fermi superfluid in the crossover region from 
a knowledge of the same of a small system through a 
scaling relation, cf. Eq. ([8]). The thus predicted energy 
and chemical potential of a system with 10 5 atoms from 
a knowledge of the same with 4 (10) atoms is found to 
have an error of less than 3% (2%). Actually, for small 
systems one needs to introduce finite-size effects, such as, 
higher-order gradient corrections and shell effects. The 
2 — 3% discrepancy quantifies the contribution of such ef- 
fects. Although we used a DF formulation in the present 
study, due to the dominance of Pauli repulsion and exis- 
tence of robust scalings, we do not believe our conclusion 
to be so peculiar as to have no general validity. 
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